SUBROUTINE DG_GRID
!-----------------------------------------------------------------------
!
! Definition von DG:  Gittergleichung
!

      use config,  only : np, BCflag
      use global,  only : tst, zz
      use primvar, only : DG, X, XZ, XA, MR, MD, ME, MUr, &
                          DG_i, DG_im1, DG_ip1, DG_im2, DG_ip2
      use physco,  only : z1, z2, z12, pi, sqrt2pi
      use gridvar, only : beta, GW_MD, GW_ME, GW_OPA, tau, &
                          n, nA, n_tmp, Res
      use zvar,    only : H_p, dH_pdE, dH_pdR
      use matvar,  only : OPAros, dOPArosdE, dOPArosdD


      implicit none

      integer :: i


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

   ! innere Pseudozellen & innere Randbedingung: i=1, i=2, i=3
      DG(MR,MR,DG_i,1)      =  z1
      DG(MR,MR,DG_i,2)      =  z1
      DG(MR,MR,DG_i,3)      =  z1

   ! aeussere Pseudozellen & aeussere Randbedingung: i=np, i=np-1, i=np-2

      DG(MR,MR,DG_i,np)     =  z1
      DG(MR,MR,DG_i,np-1)   =  z1
      DG(MR,MR,DG_i,np-2)   =  z1

   ! innerster Gitterpunkt IM Gitter: i=4 - konstante Gitterkonzentration

DG(MR,MR,DG_im2,4)=beta*(tau(4-1) / tst + z1 )*(((X(MR,4-1) + X(MR,4-2))*(z1  - BCflag(4-2))*z12 ) / &
(-(X(MR,4-2)) + X(MR,4-1)*(BCflag(4-2) + z1 ))**2 + ((z1  - BCflag(4-2))*z12 ) / (-(X(MR,4-2)) + X(MR,4-1)*(BCflag(4-2) + z1 )))
DG(MR,MR,DG_im1,4)=-(beta*(tau(4) / tst + z1 )*(((X(MR,4) + X(MR,4-1))*(z1  - BCflag(4-1))*z12 ) / (-(X(MR,4-1)) &
+ X(MR,4)*(BCflag(4-1) + z1 ))**2 + ((z1  - BCflag(4-1))*z12 ) / (-(X(MR,4-1)) + X(MR,4)*(BCflag(4-1) + z1 )))) - (tau(4-1) / tst &
+ z1 )*(-(beta*(-(((X(MR,4-1) + X(MR,4-2))*(BCflag(4-2) + z1 )*(z1  - BCflag(4-2))*z12 ) / (-(X(MR,4-2)) + X(MR,4-1)*(BCflag(4-2) &
+ z1 ))**2) + ((z1  - BCflag(4-2))*z12 ) / (-(X(MR,4-2)) + X(MR,4-1)*(BCflag(4-2) + z1 )))) + ((X(MR,4) + X(MR,4-1))*z12 *(z1  + &
beta*(-(BCflag(4)) - BCflag(4-2) + z2 ))) / (-(X(MR,4-1)) + X(MR,4)*(BCflag(4-1) + z1 ))**2 + (z12 *(z1  + beta*(-(BCflag(4)) - &
BCflag(4-2) + z2 ))) / (-(X(MR,4-1)) + X(MR,4)*(BCflag(4-1) + z1 )))
DG(MR,MR,DG_i,4)=-((tau(4-1) / tst + z1 )*(-(beta*(((X(MR,4) + X(MR,4+1))*(z1  - BCflag(4))*z12 ) / (-(X(MR,4)) &
+ X(MR,4+1)*(BCflag(4) + z1 ))**2 + ((z1  - BCflag(4))*z12 ) / (-(X(MR,4)) + X(MR,4+1)*(BCflag(4) + z1 )))) - ((X(MR,4) + &
X(MR,4-1))*(BCflag(4-1) + z1 )*z12 *(z1  + beta*(-(BCflag(4)) - BCflag(4-2) + z2 ))) / (-(X(MR,4-1)) + X(MR,4)*(BCflag(4-1) + z1 &
))**2 + (z12 *(z1  + beta*(-(BCflag(4)) - BCflag(4-2) + z2 ))) / (-(X(MR,4-1)) + X(MR,4)*(BCflag(4-1) + z1 )))) + (tau(4) / tst + &
z1 )*(-(beta*(-(((X(MR,4) + X(MR,4-1))*(BCflag(4-1) + z1 )*(z1  - BCflag(4-1))*z12 ) / (-(X(MR,4-1)) + X(MR,4)*(BCflag(4-1) + z1 &
))**2) + ((z1  - BCflag(4-1))*z12 ) / (-(X(MR,4-1)) + X(MR,4)*(BCflag(4-1) + z1 )))) + ((X(MR,4) + X(MR,4+1))*z12 *(z1  + &
beta*(-(BCflag(4-1)) - BCflag(4+1) + z2 ))) / (-(X(MR,4)) + X(MR,4+1)*(BCflag(4) + z1 ))**2 + (z12 *(z1  + beta*(-(BCflag(4-1)) - &
BCflag(4+1) + z2 ))) / (-(X(MR,4)) + X(MR,4+1)*(BCflag(4) + z1 )))
DG(MR,MR,DG_ip1,4)=beta*(tau(4-1) / tst + z1 )*(-(((X(MR,4) + X(MR,4+1))*(BCflag(4) + z1 )*(z1  - BCflag(4))*z12 &
) / (-(X(MR,4)) + X(MR,4+1)*(BCflag(4) + z1 ))**2) + ((z1  - BCflag(4))*z12 ) / (-(X(MR,4)) + X(MR,4+1)*(BCflag(4) + z1 ))) + &
(tau(4) / tst + z1 )*(-(beta*(((X(MR,4+1) + X(MR,4+2))*(z1  - BCflag(4+1))*z12 ) / (-(X(MR,4+1)) + X(MR,4+2)*(BCflag(4+1) + z1 &
))**2 + ((z1  - BCflag(4+1))*z12 ) / (-(X(MR,4+1)) + X(MR,4+2)*(BCflag(4+1) + z1 )))) - ((X(MR,4) + X(MR,4+1))*(BCflag(4) + z1 &
)*z12 *(z1  + beta*(-(BCflag(4-1)) - BCflag(4+1) + z2 ))) / (-(X(MR,4)) + X(MR,4+1)*(BCflag(4) + z1 ))**2 + (z12 *(z1  + &
beta*(-(BCflag(4-1)) - BCflag(4+1) + z2 ))) / (-(X(MR,4)) + X(MR,4+1)*(BCflag(4) + z1 )))
DG(MR,MR,DG_ip2,4)=-(beta*(tau(4) / tst + z1 )*(-(((X(MR,4+1) + X(MR,4+2))*(BCflag(4+1) + z1 )*(z1  - &
BCflag(4+1))*z12 ) / (-(X(MR,4+1)) + X(MR,4+2)*(BCflag(4+1) + z1 ))**2) + ((z1  - BCflag(4+1))*z12 ) / (-(X(MR,4+1)) + &
X(MR,4+2)*(BCflag(4+1) + z1 ))))

!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=5,np-3

DG(MR,MR,DG_im2,i)=(beta*(tau(i-1) / tst + z1 )*(((X(MR,i-1) + X(MR,i-2))*(z1  - BCflag(i-2))*z12 ) / &
(-(X(MR,i-2)) + X(MR,i-1)*(BCflag(i-2) + z1 ))**2 + ((z1  - BCflag(i-2))*z12 ) / (-(X(MR,i-2)) + X(MR,i-1)*(BCflag(i-2) + z1 )))) &
/ Sqrt(z1  + ((X(MR,i) + X(MR,i-1))**2*z12 **2*(GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + &
GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / &
OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2) + (n_tmp(i-1)*(X(MR,i) + &
X(MR,i-1))**2*z12 **3*((GW_OPA*X(MD,i-2)*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - &
OPAros(i-2))*dH_pdR(i-2)*dOPArosdD(i-2)) / (sqrt2pi*H_p(i-2)**2) + (GW_OPA*X(MD,i-2)*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / &
OPAros(i-2))*(OPAros(i-1) - OPAros(i-2))**2*dH_pdR(i-2)*dOPArosdD(i-2)) / (sqrt2pi*H_p(i-2)**2*OPAros(i-2)**2))) / ((-(X(MR,i-1)) &
+ X(MR,i)*(BCflag(i-1) + z1 ))**2*(z1  + ((X(MR,i) + X(MR,i-1))**2*z12 **2*(GW_ME*(z1  / X(ME,i-1) + z1  / &
X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 &
**2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / (-(X(MR,i-1)) + &
X(MR,i)*(BCflag(i-1) + z1 ))**2)**1.5)
DG(MR,MR,DG_im1,i)=-((beta*(tau(i) / tst + z1 )*(((X(MR,i) + X(MR,i-1))*(z1  - BCflag(i-1))*z12 ) / &
(-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2 + ((z1  - BCflag(i-1))*z12 ) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 )))) / &
Sqrt(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 **2*(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + &
GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / &
OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2)) - ((tau(i-1) / tst + z1 &
)*(-(beta*(-(((X(MR,i-1) + X(MR,i-2))*(BCflag(i-2) + z1 )*(z1  - BCflag(i-2))*z12 ) / (-(X(MR,i-2)) + X(MR,i-1)*(BCflag(i-2) + z1 &
))**2) + ((z1  - BCflag(i-2))*z12 ) / (-(X(MR,i-2)) + X(MR,i-1)*(BCflag(i-2) + z1 )))) + ((X(MR,i) + X(MR,i-1))*z12 *(z1  + &
beta*(-(BCflag(i)) - BCflag(i-2) + z2 ))) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2 + (z12 *(z1  + beta*(-(BCflag(i)) - &
BCflag(i-2) + z2 ))) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 )))) / Sqrt(z1  + ((X(MR,i) + X(MR,i-1))**2*z12 **2*(GW_ME*(z1  &
/ X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - &
X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / &
(-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2) - (n_tmp(i)*(X(MR,i) + X(MR,i+1))**2*z12 **3*((GW_OPA*X(MD,i-1)*z12 **2*z2 *(z1  &
/ OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))*dH_pdR(i-1)*dOPArosdD(i-1)) / (sqrt2pi*H_p(i-1)**2) + &
(GW_OPA*X(MD,i-1)*z12 **2*z2 *(z1  / OPAros(i) + z1  / OPAros(i-1))*(OPAros(i) - OPAros(i-1))**2*dH_pdR(i-1)*dOPArosdD(i-1)) / &
(sqrt2pi*H_p(i-1)**2*OPAros(i-1)**2))) / ((-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2*(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 &
**2*(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + GW_MD*(z1  / X(MD,i) + z1  / &
X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - &
OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2)**1.5) + (n_tmp(i-1)*z12 *(((X(MR,i) + X(MR,i-1))**2*z12 **2*z2 &
*(GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / &
X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - &
OPAros(i-2))**2)) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**3 + ((X(MR,i) + X(MR,i-1))*z12 **2*z2 *(GW_ME*(z1  / X(ME,i-1) &
+ z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - &
X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / &
(-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2 + ((X(MR,i) + X(MR,i-1))**2*z12 **2*(-((GW_OPA*X(MD,i-1)*z12 **2*z2 *(z1  / &
OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))*dH_pdR(i-1)*dOPArosdD(i-1)) / (sqrt2pi*H_p(i-1)**2)) + &
(GW_OPA*X(MD,i-1)*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / OPAros(i-2))*(OPAros(i-1) - OPAros(i-2))**2*dH_pdR(i-1)*dOPArosdD(i-1)) &
/ (sqrt2pi*H_p(i-1)**2*OPAros(i-1)**2))) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2)) / (z1  + ((X(MR,i) + &
X(MR,i-1))**2*z12 **2*(GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + &
z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - &
OPAros(i-2))**2)) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2)**1.5
DG(MR,MR,DG_i,i)=((tau(i) / tst + z1 )*(-(beta*(-(((X(MR,i) + X(MR,i-1))*(BCflag(i-1) + z1 )*(z1  - &
BCflag(i-1))*z12 ) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2) + ((z1  - BCflag(i-1))*z12 ) / (-(X(MR,i-1)) + &
X(MR,i)*(BCflag(i-1) + z1 )))) + ((X(MR,i) + X(MR,i+1))*z12 *(z1  + beta*(-(BCflag(i-1)) - BCflag(i+1) + z2 ))) / (-(X(MR,i)) + &
X(MR,i+1)*(BCflag(i) + z1 ))**2 + (z12 *(z1  + beta*(-(BCflag(i-1)) - BCflag(i+1) + z2 ))) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + &
z1 )))) / Sqrt(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 **2*(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 &
**2 + GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / &
OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2) - ((tau(i-1) / tst + z1 &
)*(-(beta*(((X(MR,i) + X(MR,i+1))*(z1  - BCflag(i))*z12 ) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2 + ((z1  - &
BCflag(i))*z12 ) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 )))) - ((X(MR,i) + X(MR,i-1))*(BCflag(i-1) + z1 )*z12 *(z1  + &
beta*(-(BCflag(i)) - BCflag(i-2) + z2 ))) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2 + (z12 *(z1  + beta*(-(BCflag(i)) - &
BCflag(i-2) + z2 ))) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 )))) / Sqrt(z1  + ((X(MR,i) + X(MR,i-1))**2*z12 **2*(GW_ME*(z1  &
/ X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - &
X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / &
(-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2) + (n_tmp(i-1)*z12 *(-(((X(MR,i) + X(MR,i-1))**2*(BCflag(i-1) + z1 )*z12 **2*z2 &
*(GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / &
X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - &
OPAros(i-2))**2)) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**3) + ((X(MR,i) + X(MR,i-1))*z12 **2*z2 *(GW_ME*(z1  / X(ME,i-1) &
+ z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - &
X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / &
(-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2)) / (z1  + ((X(MR,i) + X(MR,i-1))**2*z12 **2*(GW_ME*(z1  / X(ME,i-1) + z1  / &
X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 &
**2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / (-(X(MR,i-1)) + &
X(MR,i)*(BCflag(i-1) + z1 ))**2)**1.5 - (n_tmp(i)*z12 *(((X(MR,i) + X(MR,i+1))**2*z12 **2*z2 *(GW_ME*(z1  / X(ME,i) + z1  / &
X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + &
GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 &
))**3 + ((X(MR,i) + X(MR,i+1))*z12 **2*z2 *(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + &
GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / &
OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2 + ((X(MR,i) + X(MR,i+1))**2*z12 &
**2*(-((GW_OPA*X(MD,i)*z12 **2*z2 *(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))*dH_pdR(i)*dOPArosdD(i)) / &
(sqrt2pi*H_p(i)**2)) + (GW_OPA*X(MD,i)*z12 **2*z2 *(z1  / OPAros(i) + z1  / OPAros(i-1))*(OPAros(i) - &
OPAros(i-1))**2*dH_pdR(i)*dOPArosdD(i)) / (sqrt2pi*H_p(i)**2*OPAros(i)**2))) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2)) / &
(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 **2*(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + GW_MD*(z1 &
 / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / &
OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2)**1.5
DG(MR,MR,DG_ip1,i)=((tau(i) / tst + z1 )*(-(beta*(((X(MR,i+1) + X(MR,i+2))*(z1  - BCflag(i+1))*z12 ) / &
(-(X(MR,i+1)) + X(MR,i+2)*(BCflag(i+1) + z1 ))**2 + ((z1  - BCflag(i+1))*z12 ) / (-(X(MR,i+1)) + X(MR,i+2)*(BCflag(i+1) + z1 )))) &
- ((X(MR,i) + X(MR,i+1))*(BCflag(i) + z1 )*z12 *(z1  + beta*(-(BCflag(i-1)) - BCflag(i+1) + z2 ))) / (-(X(MR,i)) + &
X(MR,i+1)*(BCflag(i) + z1 ))**2 + (z12 *(z1  + beta*(-(BCflag(i-1)) - BCflag(i+1) + z2 ))) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + &
z1 )))) / Sqrt(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 **2*(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 &
**2 + GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / &
OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2) - (n_tmp(i)*z12 *(-(((X(MR,i) + &
X(MR,i+1))**2*(BCflag(i) + z1 )*z12 **2*z2 *(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + &
GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / &
OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**3) + ((X(MR,i) + X(MR,i+1))*z12 &
**2*z2 *(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + GW_MD*(z1  / X(MD,i) + z1  / &
X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - &
OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2)) / (z1  + ((X(MR,i) + X(MR,i+1))**2*z12 **2*(GW_ME*(z1  / &
X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - &
X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + &
X(MR,i+1)*(BCflag(i) + z1 ))**2)**1.5 + (beta*(tau(i-1) / tst + z1 )*(-(((X(MR,i) + X(MR,i+1))*(BCflag(i) + z1 )*(z1  - &
BCflag(i))*z12 ) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2) + ((z1  - BCflag(i))*z12 ) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) &
+ z1 )))) / Sqrt(z1  + ((X(MR,i) + X(MR,i-1))**2*z12 **2*(GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - &
X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / &
OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2)
DG(MR,MR,DG_ip2,i)=-((beta*(tau(i) / tst + z1 )*(-(((X(MR,i+1) + X(MR,i+2))*(BCflag(i+1) + z1 )*(z1  - &
BCflag(i+1))*z12 ) / (-(X(MR,i+1)) + X(MR,i+2)*(BCflag(i+1) + z1 ))**2) + ((z1  - BCflag(i+1))*z12 ) / (-(X(MR,i+1)) + &
X(MR,i+2)*(BCflag(i+1) + z1 )))) / Sqrt(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 **2*(GW_ME*(z1  / X(ME,i) + z1  / &
X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + &
GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 &
))**2))
DG(MR,MD,DG_im2,i)=(n_tmp(i-1)*(X(MR,i) + X(MR,i-1))**2*z12 **3*(-(GW_MD*(z1  / X(MD,i-1) + z1  / &
X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))*z12 **2*z2 ) - (GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))*(X(MD,i-1) - X(MD,i-2))**2*z12 &
**2*z2 ) / X(MD,i-2)**2 - (GW_OPA*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - &
OPAros(i-2))*dOPArosdD(i-2)) / (sqrt2pi*H_p(i-2)) - (GW_OPA*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / OPAros(i-2))*(OPAros(i-1) - &
OPAros(i-2))**2*dOPArosdD(i-2)) / (sqrt2pi*H_p(i-2)*OPAros(i-2)**2))) / ((-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2*(z1  + &
((X(MR,i) + X(MR,i-1))**2*z12 **2*(GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / &
X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / &
OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2)**1.5)
DG(MR,MD,DG_im1,i)=-((n_tmp(i)*(X(MR,i) + X(MR,i+1))**2*z12 **3*(-(GW_MD*(z1  / X(MD,i) + z1  / &
X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))*z12 **2*z2 ) - (GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))*(X(MD,i) - X(MD,i-1))**2*z12 **2*z2 ) &
/ X(MD,i-1)**2 - (GW_OPA*z12 **2*z2 *(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))*dOPArosdD(i-1)) / &
(sqrt2pi*H_p(i-1)) - (GW_OPA*z12 **2*z2 *(z1  / OPAros(i) + z1  / OPAros(i-1))*(OPAros(i) - OPAros(i-1))**2*dOPArosdD(i-1)) / &
(sqrt2pi*H_p(i-1)*OPAros(i-1)**2))) / ((-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2*(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 &
**2*(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + GW_MD*(z1  / X(MD,i) + z1  / &
X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - &
OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2)**1.5)) + (n_tmp(i-1)*(X(MR,i) + X(MR,i-1))**2*z12 &
**3*(GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))*z12 **2*z2  - (GW_MD*(z1  / X(MD,i-1) + z1  / &
X(MD,i-2))*(X(MD,i-1) - X(MD,i-2))**2*z12 **2*z2 ) / X(MD,i-1)**2 + (GW_OPA*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / &
OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))*dOPArosdD(i-1)) / (sqrt2pi*H_p(i-1)) - (GW_OPA*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / &
OPAros(i-2))*(OPAros(i-1) - OPAros(i-2))**2*dOPArosdD(i-1)) / (sqrt2pi*H_p(i-1)*OPAros(i-1)**2))) / ((-(X(MR,i-1)) + &
X(MR,i)*(BCflag(i-1) + z1 ))**2*(z1  + ((X(MR,i) + X(MR,i-1))**2*z12 **2*(GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) &
- X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / &
OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))**2)) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2)**1.5)
DG(MR,MD,DG_i,i)=-((n_tmp(i)*(X(MR,i) + X(MR,i+1))**2*z12 **3*(GW_MD*(z1  / X(MD,i) + z1  / &
X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))*z12 **2*z2  - (GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))*(X(MD,i) - X(MD,i-1))**2*z12 **2*z2 ) &
/ X(MD,i)**2 + (GW_OPA*z12 **2*z2 *(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))*dOPArosdD(i)) / &
(sqrt2pi*H_p(i)) - (GW_OPA*z12 **2*z2 *(z1  / OPAros(i) + z1  / OPAros(i-1))*(OPAros(i) - OPAros(i-1))**2*dOPArosdD(i)) / &
(sqrt2pi*H_p(i)*OPAros(i)**2))) / ((-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2*(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 &
**2*(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + GW_MD*(z1  / X(MD,i) + z1  / &
X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - &
OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2)**1.5))
DG(MR,ME,DG_im2,i)=(n_tmp(i-1)*(X(MR,i) + X(MR,i-1))**2*z12 **3*(-(GW_ME*(z1  / X(ME,i-1) + z1  / &
X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))*z12 **2*z2 ) - (GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))*(X(ME,i-1) - X(ME,i-2))**2*z12 &
**2*z2 ) / X(ME,i-2)**2 - (GW_OPA*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / OPAros(i-2))*(OPAros(i-1) - &
OPAros(i-2))**2*(-((X(MD,i-2)*dOPArosdD(i-2)*dH_pdE(i-2)) / (sqrt2pi*H_p(i-2)**2)) + dOPArosdE(i-2))) / OPAros(i-2)**2 + &
GW_OPA*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - OPAros(i-2))*((X(MD,i-2)*dOPArosdD(i-2)*dH_pdE(i-2)) &
/ (sqrt2pi*H_p(i-2)**2) - dOPArosdE(i-2)))) / ((-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2*(z1  + ((X(MR,i) + &
X(MR,i-1))**2*z12 **2*(GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) + &
z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - &
OPAros(i-2))**2)) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2)**1.5)
DG(MR,ME,DG_im1,i)=(n_tmp(i-1)*(X(MR,i) + X(MR,i-1))**2*z12 **3*(GW_ME*(z1  / X(ME,i-1) + z1  / &
X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))*z12 **2*z2  - (GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))*(X(ME,i-1) - X(ME,i-2))**2*z12 &
**2*z2 ) / X(ME,i-1)**2 + GW_OPA*z12 **2*z2 *(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) - &
OPAros(i-2))*(-((X(MD,i-1)*dOPArosdD(i-1)*dH_pdE(i-1)) / (sqrt2pi*H_p(i-1)**2)) + dOPArosdE(i-1)) - (GW_OPA*z12 **2*z2 *(z1  / &
OPAros(i-1) + z1  / OPAros(i-2))*(OPAros(i-1) - OPAros(i-2))**2*(-((X(MD,i-1)*dOPArosdD(i-1)*dH_pdE(i-1)) / &
(sqrt2pi*H_p(i-1)**2)) + dOPArosdE(i-1))) / OPAros(i-1)**2)) / ((-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2*(z1  + ((X(MR,i) &
+ X(MR,i-1))**2*z12 **2*(GW_ME*(z1  / X(ME,i-1) + z1  / X(ME,i-2))**2*(X(ME,i-1) - X(ME,i-2))**2*z12 **2 + GW_MD*(z1  / X(MD,i-1) &
+ z1  / X(MD,i-2))**2*(X(MD,i-1) - X(MD,i-2))**2*z12 **2 + GW_OPA*z12 **2*(z1  / OPAros(i-1) + z1  / OPAros(i-2))**2*(OPAros(i-1) &
- OPAros(i-2))**2)) / (-(X(MR,i-1)) + X(MR,i)*(BCflag(i-1) + z1 ))**2)**1.5) - (n_tmp(i)*(X(MR,i) + X(MR,i+1))**2*z12 &
**3*(-(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))*z12 **2*z2 ) - (GW_ME*(z1  / X(ME,i) + z1  / &
X(ME,i-1))*(X(ME,i) - X(ME,i-1))**2*z12 **2*z2 ) / X(ME,i-1)**2 - (GW_OPA*z12 **2*z2 *(z1  / OPAros(i) + z1  / &
OPAros(i-1))*(OPAros(i) - OPAros(i-1))**2*(-((X(MD,i-1)*dOPArosdD(i-1)*dH_pdE(i-1)) / (sqrt2pi*H_p(i-1)**2)) + dOPArosdE(i-1))) / &
OPAros(i-1)**2 + GW_OPA*z12 **2*z2 *(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - &
OPAros(i-1))*((X(MD,i-1)*dOPArosdD(i-1)*dH_pdE(i-1)) / (sqrt2pi*H_p(i-1)**2) - dOPArosdE(i-1)))) / ((-(X(MR,i)) + &
X(MR,i+1)*(BCflag(i) + z1 ))**2*(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 **2*(GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))**2*(X(ME,i) - &
X(ME,i-1))**2*z12 **2 + GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 **2 + GW_OPA*z12 **2*(z1  / &
OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2)**1.5)
DG(MR,ME,DG_i,i)=-((n_tmp(i)*(X(MR,i) + X(MR,i+1))**2*z12 **3*(GW_ME*(z1  / X(ME,i) + z1  / &
X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))*z12 **2*z2  - (GW_ME*(z1  / X(ME,i) + z1  / X(ME,i-1))*(X(ME,i) - X(ME,i-1))**2*z12 **2*z2 ) &
/ X(ME,i)**2 + GW_OPA*z12 **2*z2 *(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - &
OPAros(i-1))*(-((X(MD,i)*dOPArosdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dOPArosdE(i)) - (GW_OPA*z12 **2*z2 *(z1  / OPAros(i) + &
z1  / OPAros(i-1))*(OPAros(i) - OPAros(i-1))**2*(-((X(MD,i)*dOPArosdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dOPArosdE(i))) / &
OPAros(i)**2)) / ((-(X(MR,i)) + X(MR,i+1)*(BCflag(i) + z1 ))**2*(z1  + ((X(MR,i) + X(MR,i+1))**2*z12 **2*(GW_ME*(z1  / X(ME,i) + &
z1  / X(ME,i-1))**2*(X(ME,i) - X(ME,i-1))**2*z12 **2 + GW_MD*(z1  / X(MD,i) + z1  / X(MD,i-1))**2*(X(MD,i) - X(MD,i-1))**2*z12 &
**2 + GW_OPA*z12 **2*(z1  / OPAros(i) + z1  / OPAros(i-1))**2*(OPAros(i) - OPAros(i-1))**2)) / (-(X(MR,i)) + X(MR,i+1)*(BCflag(i) &
+ z1 ))**2)**1.5))

      end do


END SUBROUTINE DG_GRID

